Multifractal distribution of spike intervals for two neurons with unreliable synapses 
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Two neurons coupled by unreliable synapses are modeled by leaky integrate-and-fire neurons and 
stochastic on-off synapses. The dynamics is mapped to an iterated function system. Numerical 
calculations yield a multifractal distribution of interspike intervals. The Haussdorf, entropy and 
correlation dimensions are calculated as a function of synaptic strength and transmission probability. 
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Neurons communicate via synaptic contacts. When a 
neuron fires it sends an electric pulse (spike) along its 
axon. This spike activates biochemical processes in the 
vicinity of a synaptic contact which change the electrical 
membrane potential at the neighboring neuron. How- 
ever, experiments on synaptic contacts show that this 
complex biophysical and biochemical process is not de- 
terministic. Any incoming electrical pulse activates the 
synapse with some probability, only. In the cortex, trans- 
mission probabilities between 10% and 90% are reported 
P, 0| ■ Although model calculations show that stochas- 
tic synapses can transmit information Q it is still an 
unsolved mystery how a neural network with unreliable 
synapses is able to perform reliable computations. 

A quantitative measure of the activity of neurons is 
the distribution of interspike intervals. Typically, one 
observes broad distributions which may be described by 
a simple mathematical approach: Each neuron is mod- 
eled by a stochastic process which is driven by random 
uncorrelated synaptic inputs. Hence, usually the effect of 
unreliable synapses is modeled by external uncorrelated 
noise 0,0. 

In this paper we investigate the dynamics of two neu- 
rons coupled by unreliable synapses. The synapses are 
explicitly modeled by a Bernoulli process: Any synapse 
transmits the spike with some probability p which is inde- 
pendent of the state of the system. Our approach allows 
to calculate the spike intervals from an iterated-function- 
system (IFS). Our main result is a multifractal distribu- 
tion of interspike intervals. The Haussdorf, entropy and 
correlation dimensions are calculated as a function of the 
synaptic strength and the probability of synaptic trans- 
mission. We find a transition between connected and 
multifractal support of the distribution of spike intervals. 

In fact, fractal time series of neural spikes which are ob- 
served in many different biological systems have been re- 
lated to quantal neurotransmitter release . Our model 
shows that even a simple on-off synapse leads to fractal 
structures of the neural activity. However, our simple 
model makes predictions for the distribution of spike in- 
tervals of two coupled neurons but it does not explain 
the nature of fractal time series. 

The two neurons are modeled by a leaky integrate- 
and-fire mechanism working above threshold. In a more 
general framework, our model is a system of two identical 
pulse-coupled oscillators 0. Without synaptic contacts 



the neurons are deterministic and oscillate periodically, 
one obtains two intervals between the firing times of the 
two neurons. With reliable inhibitory synaptic contact, 
and without any delay of the synaptic transmission, the 
two neurons relax into a state of anti-phase oscillations 
with a single spike interval. With unreliable synapses, 
however, the system has a broad distribution of spike 
intervals which becomes multifractal in some range of 
the model parameters. 

Each neuron is described by the following differen- 
tial equation for the time-dependent membrane potential 
V(t): 



(i) 



As soon as the potential crosses a threshold value 9 it 
is reset to a value V r < 8. In addition it fires, i.e. it 
sends a spike to its neighbor which is transmitted with 
a probability p. If a spike is transmitted it reduces the 
potential of the receiving neuron by an amount J. For 
simplicity, we consider only inhibitory synapses to avoid 
an introduction of a refractory time. However, we believe 
that our main results do not depend on the details of the 
model. 

The neurons are working above threshold, 9 < fi, oth- 
erwise they would not fire at all. Hence the parameter 
fi controls the effect of any mechanism which forces the 
neurons to fire. Without synaptic couplings each neuron 
fires periodically with the period 



T = t In 



fl-6 



(2) 



Without loss of generality we set V r — 0, /i = 1 and r = 1, 
and in the following we use the parameter 8 = 0.95 which 
gives a period of T ~ 2.996 r. 

Figure Q shows the potential of the two neurons for a 
typical situation. At time t\ the neuron A fires and the 
spike is transmitted to neuron B resulting in a decrease 
of the potential by an amount J. The next firing event 
occurs at time t 2 . The time interval between firing events 
is denoted by A. Using the analytic solution of the differ- 
ential equation Q one obtains an iteration of the spike 
intervals A. For the quantity x = exp(— A) the iteration 
has the form 



x' = fi(x), ie {1,2,3,4,5} 



(3) 
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FIG. 1: Membrane potential of the two neurons. The spike 
of neuron A at time ti is transmitted to neuron B. 



where the five functions fi are selected according to the 
transmission probability p and the previous value of x. 
For the situation of Fig. I, which occurs with probability 
p (transmission), one finds 



1 - e 
x + J 



A(*) 



(4) 



With probability 1 — p (no transmission) the sum A + 
A' = T is identical to the period of unperturbed oscilla- 
tions which gives 



I - 



(5) 



Hence, two simple functions are iterated according to 
probability p of synaptic transmission. The situation be- 
comes slightly more complicated when neuron A over- 
takes neuron B, i.e. when one neuron fires twice before 
the other one is firing again. This occurs when the poten- 
tial Vs(ii+) becomes negative after neuron A has fired, 
that is when x > 1 — J. In this case one has A' = T or 



a/ = 1-0:= h(x) 



(6) 



But now A" depends on A and one finds with probabil- 
ity p 



1 



h{x) 



and with probability 1 — p 



1 



x + J 



fs(x) 



(7) 



(8) 



If the synaptic pulse J is larger than 9/(2 — 8) the same 
neuron can even fire more than twice in a row, but we do 
not consider such large unphysiological values of J. 

In summary, only five simple functions are iterated to 
calculate the distribution of spike intervals A. It is well 




(b) J = 0.25 



FIG. 2: Histogram of the spike intervals for the transmission 
probability p = 0.5 and the strength of the synaptic pulse 
J = 0.1 (a) and J = 0.25 (b). 



known that such a system (IFS) may lead to a fractal 
structure of the set of iterated values [8j . In our numeri- 
cal simulations of equations (4) to (8) we have generated 
about 10 11 spike intervals for each set of parameters. Fig- 
urc[21shows two histograms of the spike intervals for small 
and large values of J. Obviously, the distribution of spike 
intervals has a complex structure which we quantify by 
the Renyi dimensions [2j 



p- 



i r 
— m &* 

i=l 



(9) 



Here e is the size of the boxes of the histogram and 
is the normalized number of data points in the box i. 
The sum runs over all nonempty boxes. For (3=1, the 
entropy 1(1) = X)I=i-Pi m -Pi ^ s calculated. 

We consider three Renyi dimensions: The covering or 
box dimension D(0) which is usually identical to the 
Haussdorf dimension, the entropy dimension D(l) and 
the correlation dimension D(2). Figure shows that a 
plot of I ((3) versus lne yields a straight line over several 
orders of magnitude, hence the corresponding dimension 
can reliably be estimated from the slope of this line. In 
addition, we checked our results for the correlation di- 
mension by applying the software package TISEAN to 
our data 
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FIG. 4: Renyi dimensions (a) as a function of the strength J 
of the synaptic pulse (for p = 0.5) and (b) as a function of 
the transmission probability p (for J = 0.25) 



FIG. 5: The phases </> of the neurons are iterated by the two 
functions F\ (bottom) and F2 (top) shown by the dashed lines. 
The openings of the two functions show the empty intervals 
in the distribution of iterated phases. 



observes D(0) ~ D(l) ~ D(2) ~ 1. For large values 
of J the three dimensions are different, which means 
that the distribution of spike intervals is multifractal [j| . 
While the covering dimension D(0), i.e. the structure 
of the support of the distribution, is insensitive to the 
value of p, the entropy as well as the correlation dimen- 
sion decrease to the value zero in the deterministic limit 
p — ► 1. In fact, for p = 1, the distribution of spike inter- 
vals is a delta-peak at the fixed point of f\ which gives 
A = — ln(— J + \/4 + J 2 — 40)/2. Surprisingly, even for 
p < 1 the distribution has its maximum at this value, a 
sharp peak, as can be seen from Fig. [5] 

The results of Fig. 0fa) do not rule out a sharp tran- 
sition between a smooth and multifractal distribution of 
spike intervals. In fact, for the covering dimension -D(O), 
the transition point can be found analytically. It is conve- 
nient to transform Eq. 0) to d(j)/dt = 1 where the phase 
4> is defined as 



(f)(V) = -ln(l - V) 



(10) 



The results for the three different Renyi dimensions 
are shown in Fig. Of course, our results obey the 
exact relations D(2) < D{1) < D(0). With increas- 
ing coupling strength J and transmission probability p 
the three dimensions decrease. For small values of J 
the distribution of spike intervals is smooth, hence one Fi(4>) 



Now we consider the phase which one neuron occupies 
after the other one has fired. After the neuron A has fired 
it has the phase = 0, whereas the other neuron B has 
a nonzero phase (pi. If <pi is positive it will be neuron B 
which fires next, namely after the time T — fa. However, 
if fa is negative then neuron A will fire again after the 
time T. Regardless of which neuron fires, in both cases 
we record the phase fa+i of the neuron which has not 
fired. Given a phase fa, the next phase fa+i results by 
applying one of two mappings depending on whether a 
spike has been transmitted at time tj+i or not. These 
two mappings fa+i — F\{fa) and fa+i = F2(fa) which 
describe the transformation of phases are as follows (see 
Fig. : 



ln[exp(|0| — T) + J] (transmission) (11) 



4 



F 2 {<t>) = T-\4>\ (no transmission) (12) 

The function F 2 just flips the lower interval [0, T/2] to the 
upper one [T/2, T]. The function F\ maps the complete 
interval [0, T] to the interval [- ln(l + J), - ln(l - 6 + J)} . 
If the maximum of Fi is smaller than T/2, then there 
exists an interval in the vicinity of T/2 which cannot be 
reached from outside. In Fig. [5] this interval is indicated 
by the small square in the center of the figure. This 
interval in the center is either flipped by F2 or mapped 
to an interval outside of it by F\ . This means that finally 
any point inside the square will leave it. In addition, no 
other point can enter this interval. Hence the distribution 
of phases has an opening on this interval. By consecutive 
iterations of Ti and F 2 this opening is distributed on the 
complete range of phases, as depicted in Fig. [5] by the 
openings in the functions F\ and F 2 . This indicates - but 
does not prove it - that the support of the distribution of 
spike intervals has a fractal structure, leading to D(Q) < 
1. By these arguments the support of the distribution 
has a fractal structure if the maximum of Ti is smaller 
than T/2 which gives a critical point 

J* = Vl - - (1 - 6) (13) 

For J < J* the distribution fills the complete range of 
(j) values, while for large values of J the distribution has 
empty intervals. Indeed, this value is consistent with 
the data of Fig. 0Ja) where the covering dimension -D(O) 
deviates from the value D(0) = 1 at about J* = 0.1736. 



Note, however, that even below J» the distribution is 
multifractal because the values of D(l) and D{2) are still 
smaller than one. We do not know whether there is a 
sharp transition to a smooth structure for small J values 
or whether the fractal dimensions D(l) and D(2) just 
come very close to the value one. The data of Fig. 01 do 
not allow to distinguish between these two possibilities. 

Our system of two identical pulse-coupled oscillators 
with random on-off synapses is very simplified model of 
two coupled neurons. For instance, synaptic transmis- 
sion may be multi- valued [TH Il2j | and time-delayed |13| , 
and a much better model would include the dynamics of 
ion channels |l4j| . However, in any model a random un- 
corrected process which opens and closes synaptic trans- 
mission always yields an iterated function system which 
produces fractal distributions of spike intervals depend- 
ing on the model parameters. Up to now, a fractal struc- 
ture of spike intervals has not yet been observed. But, 
to our knowlege, experiments on two interacting neurons 
under controlled conditions have not yet been reported, 
either. Our model makes predictions for such an experi- 
ment which may help to clarify the stochastic nature of 
synaptic transmission. 
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